# Data were downloaded from http://diposit.ub.edu/dspace/handle/2445/173199
saint = read.csv("Ivermectin_study_data.csv", header=T)
# 1: ivermectin ; 2: placebo - we set placebo to be 1
saint$treat[saint$treat==2]=0
# take mean value of N and E PCR
saint = dplyr::filter(saint, !(is.na(ct_n) & is.na(ct_e)))
saint$ct = rowMeans(saint[, c('ct_n','ct_e')], na.rm = T)
ind_fit = saint$day<=21

iv_dat1 = data.frame(t = saint$day[ind_fit],
                     Trt = as.factor(saint$treat)[ind_fit],
                     id = saint$studyno[ind_fit],
                     y = saint$ct[ind_fit])

run the model

CT_threshold = 32
data_stan = iv_dat1[iv_dat1$t<10, ]
ind_obs = data_stan$y<CT_threshold
ind_cens = !ind_obs
data_stan$id = as.numeric(as.factor(data_stan$id))
data_stan$Trt = as.numeric(as.factor(data_stan$Trt))-1
data_1 = list(Nobs = sum(ind_obs),
              Ncens = sum(ind_cens),
              J = length(unique(data_stan$id)),
              id_obs = data_stan$id[ind_obs],
              t_obs = (data_stan$t-1)[ind_obs],
              Trt_obs = data_stan$Trt[ind_obs],
              y_obs = data_stan$y[ind_obs],
              id_cens = data_stan$id[ind_cens],
              t_cens = (data_stan$t-1)[ind_cens],
              Trt_cens = data_stan$Trt[ind_cens],
              y_cens = data_stan$y[ind_cens],
              CT_threshold=CT_threshold,
              prior_Trt_sd=.1,
              prior_CT_intercept = 20,
              prior_sigma_sd = 1,
              prior_mean_slope = 2,
              prior_tau_intercept_sd = 2,
              prior_tau_slope_sd = 1,
              prior_tau_intercept_mean = 3,
              prior_tau_slope_mean = .5)
interim_fit1 = sampling(object = stan_lin_model,
                        data = data_1,
                        chains=4, iter=10^5,thin=10)
## 
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:     1 / 100000 [  0%]  (Warmup)
## Chain 1: Iteration: 10000 / 100000 [ 10%]  (Warmup)
## Chain 1: Iteration: 20000 / 100000 [ 20%]  (Warmup)
## Chain 1: Iteration: 30000 / 100000 [ 30%]  (Warmup)
## Chain 1: Iteration: 40000 / 100000 [ 40%]  (Warmup)
## Chain 1: Iteration: 50000 / 100000 [ 50%]  (Warmup)
## Chain 1: Iteration: 50001 / 100000 [ 50%]  (Sampling)
## Chain 1: Iteration: 60000 / 100000 [ 60%]  (Sampling)
## Chain 1: Iteration: 70000 / 100000 [ 70%]  (Sampling)
## Chain 1: Iteration: 80000 / 100000 [ 80%]  (Sampling)
## Chain 1: Iteration: 90000 / 100000 [ 90%]  (Sampling)
## Chain 1: Iteration: 100000 / 100000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 17.8136 seconds (Warm-up)
## Chain 1:                21.2163 seconds (Sampling)
## Chain 1:                39.0299 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:     1 / 100000 [  0%]  (Warmup)
## Chain 2: Iteration: 10000 / 100000 [ 10%]  (Warmup)
## Chain 2: Iteration: 20000 / 100000 [ 20%]  (Warmup)
## Chain 2: Iteration: 30000 / 100000 [ 30%]  (Warmup)
## Chain 2: Iteration: 40000 / 100000 [ 40%]  (Warmup)
## Chain 2: Iteration: 50000 / 100000 [ 50%]  (Warmup)
## Chain 2: Iteration: 50001 / 100000 [ 50%]  (Sampling)
## Chain 2: Iteration: 60000 / 100000 [ 60%]  (Sampling)
## Chain 2: Iteration: 70000 / 100000 [ 70%]  (Sampling)
## Chain 2: Iteration: 80000 / 100000 [ 80%]  (Sampling)
## Chain 2: Iteration: 90000 / 100000 [ 90%]  (Sampling)
## Chain 2: Iteration: 100000 / 100000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 18.438 seconds (Warm-up)
## Chain 2:                21.152 seconds (Sampling)
## Chain 2:                39.59 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:     1 / 100000 [  0%]  (Warmup)
## Chain 3: Iteration: 10000 / 100000 [ 10%]  (Warmup)
## Chain 3: Iteration: 20000 / 100000 [ 20%]  (Warmup)
## Chain 3: Iteration: 30000 / 100000 [ 30%]  (Warmup)
## Chain 3: Iteration: 40000 / 100000 [ 40%]  (Warmup)
## Chain 3: Iteration: 50000 / 100000 [ 50%]  (Warmup)
## Chain 3: Iteration: 50001 / 100000 [ 50%]  (Sampling)
## Chain 3: Iteration: 60000 / 100000 [ 60%]  (Sampling)
## Chain 3: Iteration: 70000 / 100000 [ 70%]  (Sampling)
## Chain 3: Iteration: 80000 / 100000 [ 80%]  (Sampling)
## Chain 3: Iteration: 90000 / 100000 [ 90%]  (Sampling)
## Chain 3: Iteration: 100000 / 100000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 17.9827 seconds (Warm-up)
## Chain 3:                21.4706 seconds (Sampling)
## Chain 3:                39.4533 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'f8a746a2461316d8be14aecf5c16c64f' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:     1 / 100000 [  0%]  (Warmup)
## Chain 4: Iteration: 10000 / 100000 [ 10%]  (Warmup)
## Chain 4: Iteration: 20000 / 100000 [ 20%]  (Warmup)
## Chain 4: Iteration: 30000 / 100000 [ 30%]  (Warmup)
## Chain 4: Iteration: 40000 / 100000 [ 40%]  (Warmup)
## Chain 4: Iteration: 50000 / 100000 [ 50%]  (Warmup)
## Chain 4: Iteration: 50001 / 100000 [ 50%]  (Sampling)
## Chain 4: Iteration: 60000 / 100000 [ 60%]  (Sampling)
## Chain 4: Iteration: 70000 / 100000 [ 70%]  (Sampling)
## Chain 4: Iteration: 80000 / 100000 [ 80%]  (Sampling)
## Chain 4: Iteration: 90000 / 100000 [ 90%]  (Sampling)
## Chain 4: Iteration: 100000 / 100000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 17.9634 seconds (Warm-up)
## Chain 4:                22.2008 seconds (Sampling)
## Chain 4:                40.1642 seconds (Total)
## Chain 4:
out=extract(interim_fit1)
par(las=1, bty='n', family='serif', 
    mfrow=c(1,2), cex.lab=1.5, cex.axis=1.5)
xs = jitter(iv_dat1$t, amount = .3)-1
plot(xs, iv_dat1$y, pch=20, xlab = 'Time from enrollment',
     ylab='Cycle threshold', xaxt='n', ylim=rev(range(iv_dat1$y)))
axis(1, at=c(0,3,6,13,20))
abline(h=32, lty=2, lwd=3, v=7, col='red')
for(id in unique(iv_dat1$id)){
  lines(xs[iv_dat1$id==id], iv_dat1$y[iv_dat1$id==id],lty=2)
}
mtext(text = 'A)', side = 3, cex=1.5, adj = 0)

hist(100*out$beta_Trt, xlab='Treatment effect (%)',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(-300:300, dnorm(c(-300:300), mean = 0, sd = 100*data_1$prior_Trt_sd), col='blue',lwd=3)
abline(v=0, lwd=3, col='red')
mtext(text = 'B)', side = 3, cex=1.5, adj = 0)

writeLines(sprintf('Probability that treatment effect>0 is %s', mean(out$beta_Trt > 0)))
## Probability that treatment effect>0 is 0.74765
par(mfrow=c(5,5), mar=c(3,3,1,0),las=1, family='serif', cex.lab=1.5, cex.axis=1.5)
for(id in 1:data_1$J){
  indo = data_1$id_obs==id
  indc = data_1$id_cens==id
  plot(data_1$t_obs[indo], 
       data_1$y_obs[indo],pch=20, xlab='', ylab='',cex=2,
       ylim= rev(range(saint$ct)), xlim = c(0,6), yaxt='n', xaxt='n')
  points(data_1$t_cens[indc], 
       data_1$y_cens[indc],pch=18,cex=2,col='blue')
  polygon(c(c(data_1$t_obs[indo],data_1$t_cens[indc]),
            rev(c(data_1$t_obs[indo],data_1$t_cens[indc]))),
          c(c(apply(out$mu_obs,2,quantile,.1)[indo],
          apply(out$mu_cens,2,quantile,.1)[indc]), 
          rev(c(apply(out$mu_obs,2,quantile,.9)[indo],
          apply(out$mu_cens,2,quantile,.9)[indc]))),
          border = NA, col = adjustcolor('grey',.3))
  lines(c(data_1$t_obs[indo],data_1$t_cens[indc]),
        c(apply(out$mu_obs,2,mean)[indo],
          apply(out$mu_cens,2,mean)[indc]),lwd=3, col=adjustcolor('grey',.5))
  points(data_1$t_obs[indo],data_1$y_obs[indo],pch=20,cex=2)
  points(data_1$t_cens[indc],data_1$y_cens[indc],pch=18,cex=2,col='blue')

  abline(h=CT_threshold, lty=2)
  if(id%%5 == 1){ axis(2, at = seq(15,40,by=5))}
  if(id>19) {axis(1, at = c(0,3,6)) }
}

Plot the posterior summaries

par(mfrow=c(2,3), cex.lab=1.5, cex.axis=1.5)
hist(out$beta_0, xlab='Population intercept',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(seq(0,30,by=.1), 
      dnorm(seq(0,30,by=.1),
            mean = data_1$prior_CT_intercept, sd = 3),
      col='blue',lwd=3)

hist(out$beta_t, xlab='Population slope',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(seq(0,10,by=.1), 
      dnorm(seq(0,10,by=.1),
            mean = data_1$prior_mean_slope, sd = 1),
      col='blue',lwd=3)

hist(out$beta_Trt, xlab='Treatment effect',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(seq(-1,1,by=.01), 
      dnorm(seq(-1,1,by=.01),
            mean = 0, sd = data_1$prior_Trt_sd),
      col='blue',lwd=3)

hist(out$tau_0, xlab='Standard deviation: random intercept',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(seq(0,10,by=.1), 
      dnorm(seq(0,10,by=.1),
            mean = data_1$prior_tau_intercept_mean, 
            sd = data_1$prior_tau_intercept_sd),
      col='blue',lwd=3)

hist(out$tau_t, xlab='Standard deviation: random slope',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(seq(-1,3,by=.01), 
      dnorm(seq(-1,3,by=.01),
            mean = data_1$prior_tau_slope_mean, 
            sd = data_1$prior_tau_slope_sd),
      col='blue',lwd=3)

hist(out$sigma, xlab='Standard deviation: residual error',
     main='', breaks = 50, yaxt='n',ylab = '', freq = F,
     border = NA, col=adjustcolor('grey',.9))
lines(seq(0,10,by=.1), 
      dnorm(seq(0,10,by=.1),
            mean = 1, sd = data_1$prior_sigma_sd),
      col='blue',lwd=3)